import pandas as pd
import numpy as np
from matplotlib import pyplot as plt
import matplotlib
import seaborn as sns
from statsmodels.tsa.arima.model import ARIMA
import statsmodels as sm
from scipy import stats

#读取文件
data = pd.read_excel('整合数据6.0.xlsx',usecols=['年末总人口(万人)','时间'])
data.index = data['时间']
data = data.drop(columns=['时间'])
data = data[0:30]
print(data.info())

#找出有缺失值的行  data[["年末总人口(万人)"]].isnull().T.any().values
def knm(df,n):
    #找出缺失值的行
    temp = df.isnull().T.any().values
    temp_df = df.copy()
    for i in range(len(temp)):
        if temp[i] == True:
            if i < n-1:                  #前n个
                 temp_df.loc[i,"年末总人口(万人)"] = df.loc[i:i+n,"年末总人口(万人)"].mean()
            elif i > len(temp) - 1 -n:     #后n个
                 temp_df.loc[i,"年末总人口(万人)"] = df.loc[i-n:i,"年末总人口(万人)"]
            else:
                 print(df.loc[i-n:i+n,"年末总人口(万人)"])
                 temp_df.loc[i,"年末总人口(万人)"] = df.loc[i-n:i+n,"年末总人口(万人)"].mean()
                 print(i-n,i+n+1)
    return temp_df
not_miss = knm(data[["年末总人口(万人)"]],2)
data["年末总人口(万人)"] = not_miss.values
print(data)

plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号
y = data.values
x = data.index
plt.plot(x,y)
plt.show()

# 数据平稳化

data["diff_1"] = data["年末总人口(万人)"].diff(1)  #一阶差分
data["diff_2"] = data["年末总人口(万人)"].diff(1)  #二阶差分

#分别画出ACF(自相关)和PACF（偏自相关）图像
from statsmodels.graphics.tsaplots import plot_acf
from statsmodels.graphics.tsaplots import plot_pacf

def diff(df,col):
    len = df[col].size
    font = {"size":15,
       "family":"fangsong"}
    matplotlib.rc("font",**font)
    matplotlib.rcParams['axes.unicode_minus']=False

    df["diff_1"] = df[col].diff(1)  #一阶差分
    df["diff_2"] = df["diff_1"].diff(1)  #二阶差分

    #平稳数据折线图
    plt.figure(figsize=(8,4))
    plt.subplot(3,1,1)
    plt.plot(df[col].values,label="源数据")
    plt.xlim(0,40)
    plt.legend()
    plt.subplot(3,1,2)
    plt.plot(df["diff_1"].values,c="darkgreen",label="一阶差分")
    plt.plot([0,40],[0,0],"--",c = "grey")
    plt.xlim(0,40)
    plt.legend()
    plt.subplot(3,1,3)
    plt.plot(df["diff_2"].values,c="tomato",label="二阶差分")
    plt.plot([0,40],[0,0],"--",c = "grey")
    plt.xlim(0,40)
    plt.legend()
    plt.show()

    #ACF PACF
    print("-"*50,"未平稳数据ACF与PACF","-"*50)
    fig = plt.figure(figsize=(12,8))
    ax1 = fig.add_subplot(211)
    fig = plot_acf(df[col], lags=len-1,ax = ax1)
    ax2 = fig.add_subplot(212)
    plot_pacf(df[col], lags=len/2-1,ax = ax2)
    plt.show()

    #一阶差分后的ACF PACF
    print("-"*50,"一阶差分数据ACF与PACF","-"*50)
    fig = plt.figure(figsize=(12,8))
    ax1 = fig.add_subplot(211)
    fig = plot_acf(df["diff_1"][1:].values, lags=len-2,ax = ax1)
    ax2 = fig.add_subplot(212)
    plot_pacf(df["diff_1"][1:], lags=len/2-2,ax = ax2)
    plt.show()

diff(data,"年末总人口(万人)")


# 数据平稳性检验 单位根检验
#未差分平稳性检测（ADF检验、单位根检验）
from  statsmodels.tsa.stattools  import  adfuller  as  ADF
print(u'原始序列的ADF检验结果为：',  ADF(data["年末总人口(万人)"]))
#返回值依次为adf、pvalue、usedlag、nobs、critical  values、icbest、regresults、resstore  p<0.05时表示稳定
#一阶差分平稳性检测（ADF检验、单位根检验）
from  statsmodels.tsa.stattools  import  adfuller  as  ADF
print(u'一阶差分序列的ADF检验结果为：',  ADF(data["diff_1"][1:]))
#返回值依次为adf、pvalue、usedlag、nobs、critical  values、icbest、regresults、resstore  p<0.05时表示稳定

#定阶
pmax  =  int(len(data["年末总人口(万人)"])/10)  #一般阶数不超过length/10
qmax  =  int(len(data["年末总人口(万人)"])/10)  #一般阶数不超过length/10
bic_matrix  =  []  #bic矩阵
for  p  in  range(pmax+1):
    tmp  =  []
    for  q  in  range(qmax+1):  #存在部分报错，所以用try来跳过报错。
        try:
            tmp.append(ARIMA(data[u"年末总人口(万人)"],order=(p,1,q)).fit().bic)
        except:
            tmp.append(None)
    bic_matrix.append(tmp)
bic_matrix  =  pd.DataFrame(bic_matrix)  #从中可以找出最小值
print(bic_matrix)

# 找出最小值位置。
aic_matrix  =  pd.DataFrame(bic_matrix)  #从中可以找出最小值
p,q  =  bic_matrix.stack().idxmin()
print(u'BIC最小的p值和q值为：%s、%s'  %(p,q))

# 模型报告
data[u'年末总人口(万人)'] = data[u'年末总人口(万人)'].astype(float)
model = ARIMA(data['年末总人口(万人)'], order=(3,1,3)).fit()  # 建立ARIMA(0, 1, 1)模型
print(model)
# 模型检验
resid=model.resid  #残差
plt.figure(figsize=(12,8))
plt.plot(resid)
plt.show()
# 残差正态性检验
print(stats.normaltest(resid))#检验序列残差是否为正态分布    pvalue=0.00028625258929196876   <  0.05  拒绝原假设 认为残差符合正太分布
print('模型报告为：\n', model.summary())
print('预测未来20年，其预测结果如下：\n', model.forecast(20))

